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In the field of thermoacoustic energy conversion, the application of numerical analysis techniques, spe¬ 
cifically computational fluid dynamics (CFD) simulations, have gained ground in recent years. Previous 
efforts have focused on single thermoacoustic couples that were subjected to the thermoacoustic effect 
through an oscillatory boundary condition. CFD simulations of an entire thermoacoustic device are com¬ 
putationally expensive and few examples exist. The present work presents an extension of a simulation of 
a whole thermoacoustic engine that also includes a refrigeration stack. Through interaction of thermally 
generated sound waves, cooling of the working gas in this stack is demonstrated. 

© 2010 Elsevier Ltd. All rights reserved. 


1. Introduction 

The basic thermodynamic cycle in thermoacoustic devices is the 
Stirling cycle, which was developed in 1816 [1]. The original 
mechanical Stirling engine utilized two pistons and a regenerative 
heat exchanger [2]. Over the course of one cycle, the working gas 
is compressed, and it then transfers thermal energy to the heat sink, 
thus maintaining a constant temperature. Next, the gas is heated at 
constant volume by the regenerator and then is heated further at the 
heat source. This heat supply occurs while the gas is allowed to 
expand and drive the power piston, again at constant temperature. 
After expansion, the gas is displaced to the heat sink, while cooling 
at constant volume by depositing heat to the regenerator, which 
stores heat between cycle segments [2]. The first application of this 
cycle as a thermoacoustic technology occurred when Ceperley rec¬ 
ognized that sound waves could replace pistons for gas compression 
and displacement [3]. 

1.1. Thermoacoustic engines 

The simplest thermoacoustic engines (TAEs) operate with a 
standing acoustic wave in a resonance tube. This wave causes pres¬ 
sure variations as well as gas displacement. In a quarter wave¬ 
length resonator comprised of a tube with one open end and one 
closed end, the pressure anti-node of the wave is located at the 
closed end, and the velocity anti-node is located at the open end. 
When a gas is subject to such a standing wave and it interacts with 
a solid wall that is subject to a (externally imposed) temperature 
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gradient, pressure disturbances can be intensified to reach strong 
amplitudes. In practice, this interaction occurs within a regenera¬ 
tive unit commonly referred to as the “stack" in standing wave de¬ 
vices. Ceramic monolith structures are often used as stacks because 
they offer a large surface-to-volume ratio and exhibit desirable 
hydraulic performance. In order for amplification to occur, this 
temperature gradient must be larger than the critical temperature 
gradient, which is related to the temperature gradient that the gas 
would experience if it were under the influence of a sound wave in 
adiabatic conditions. The expression for this critical temperature 
gradient 


was derived by Swift [4]. It depends on the operating frequency co, 
the first order pressure and velocity in the standing wave p) and uj, 
as well as the mean gas density p m and specific heat c p . 

1.2. Thermoacoustic refrigerators without moving parts 

As stated above, an imposed temperature gradient along a porous 
medium such as a ceramic monolith placed in a resonance tube can 
produce strong pressure oscillations through the thermoacoustic 
effect. Again, it is based on the transfer of heat between the working 
gas and the solid walls at two different temperature levels and pres¬ 
sures. If we allow strong sound waves to interact with a porous med¬ 
ium, this cycle is reversed; the pressure amplitude is attenuated and 
a temperature gradient can be created along this porous medium. 
This phenomenon forms the basis of thermoacoustic refrigeration. 
Simple systems utilize mechanical drivers (such as speakers) to 
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Nomenclature 


c 

specific heat 

co angular frequency 

C/d. C sf 

correlation constants 

p density 

d 

diameter 


f 

frequency 

Subscripts 

P 

pressure 

1 first order 

Q 

heat flux 

crit critical 

Re 

Reynolds number 

m mean, time averaged 

T 

temperature 

p constant pressure 

u 

velocity 

w wire 

X 

position 


X* 

position relative to left stack edge 

Superscripts 



s standing wave 

Greek symbols 


A 

wavelength 



create the pressure oscillations. However, a TAE can also be utilized 
to drive this thermoacoustic refrigerator (TAR). Such TARs are a real¬ 
ity today; however they are limited to few, specific uses, such as gas 
liquefaction. There are several examples of advanced TARs used for 
these specialized applications. Radebaugh at the National Institute 
of Standards and Technology (NIST) built a TAR with 5 W of cooling 
power at 120 K and a low temperature at no load of 90 K [5,6]. The 
main benefit is the lack of moving parts (such as seals), thus reducing 
maintenance costs. It is noteworthy that TARs can achieve these low 
temperatures in a single stage, whereas conventional vapor com¬ 
pression refrigerators (VCRs) can only achieve «230 K in a single 
stage [6]. Cryogenic cooling is not the only application for TARs. A 
practical example of this design was given by Poese with a freezer 
for ice cream storage. This small scale chiller featured an annular 
space around the regenerative unit to achieve traveling wave phas¬ 
ing [7]. The temperature gradient achievable by thermoacoustic 
cooling is limited by the critical temperature gradient. 

2. Simulation of the thermoacoustic effect 

Baseline numerical simulations of the thermoacoustic effect 
have been demonstrated several times. Early modeling efforts done 
by Cao et al. considered a reduced domain that included a single 
stack plate and sections of resonance tube excited with an oscillat¬ 
ing boundary condition (mimicking a mechanical driver) [8], This 
model design is illustrated in Fig. 1, which also shows the full de¬ 
vice. The group investigated the temperature behavior along the 
stack plate. 

These early modeling results were compared to measurements 
[9] and predictions from linear acoustic models [10]. Further ana¬ 
lytical investigations of the thermoacoustic effect were conducted 
by Ozoe et al. [11 ]. Cao’s model [8] formed the basis for several la¬ 
ter refinements. Further investigation of the temperature behavior, 
especially deviations from strictly harmonic behavior, was done by 
Worlikar et al. [12-14]. Besnoin et al. further elaborated on the 
model by including heat exchangers at either end of the stack, 
which is comprised of a single plate. This work elaborates on the 
influence of the placement of said heat exchangers relative to the 
stack plate [15,16]. The compressible Navier-Stokes equations 
were implemented into the simulation by Ishikawa et al. [17], 
Marx et al. [18-21] combined several previous research efforts 
and utilized Cao’s domain with heat exchangers. The group inves¬ 
tigated the non-linearities in the temperature profile, noting that 
the pressure shows an almost harmonic behavior with respect to 
time, and the velocity behaves harmonically. The group also iden¬ 
tified a critical stack length where non-linearities prevail, which is 
«4 times the gas’ displacement distance centered around each 


A/4 



Fig. 1 . Computational domain used by Cao et al. [8] for the investigation of the 
thermoacoustic effect. 

stack edge. In addition, the group investigated the influence of 
stack location on all parameters, including the enthalpy flux. Pic¬ 
colo et al. utilized the Cao model to investigate the influence of 
plate spacing and heat exchanger location, as well as the effect of 
the heat transfer coefficient of the stack on the heat transfer be¬ 
tween the solid and gas [22]. 

2.3. Previous whole device simulations 

Zoontjens, finally, transitioned the investigation of the single 
plate thermoacoustic couple to a computational fluid dynamics 
(CFD) simulation. Using a commercial code (Fluent), the group 
developed a mesh that included large sections of resonator and a 
single plate stack. As with all previous models, the oscillations 
were introduced by an oscillating pressure boundary condition. 
This simulation was used extensively to investigate the fluid 
behavior during the thermoacoustic cycle [23]. 

Preceding Zoontjens attempts to investigate the thermoacoustic 
effect using CFD analysis are two additional similar, but simplified, 
examples. Hanschk et al. used Fluent 4.4.4 to simulate a Rijke tube 
(which is similar to a classic thermoacoustic engine, except that 
the oscillations are created with one heated wire screen rather 
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than a stack) [24,25]. The boundary conditions were a steady 
velocity inlet and an open end, with a heated screen in the center 
[24,25], A further example of the simulation of a Rijke tube was 
performed by Entezam et al. [26]. 

Nijeholt et al. [27] used the commercial solver ANSYS CFX to 
simulate a traveling wave engine inside a Helmholtz resonator, 
mimicking the design built by Bastyr et al. [28]. Rather than mod¬ 
eling the regenerator geometrically, the group modeled this region 
with body forces to account for the viscous forces acting on the gas. 
This includes a volume porosity to account for the solid fraction in 
the regenerator region. The pressure drop across the heat exchang¬ 
ers and regenerator are accounted for by: 


A P 


c C A fm2 

Cfd+ Re 4 d w 


( 2 ) 


following the work done by Thomas et al. on the friction and heat 
transfer behavior of gases in Stirling regenerators [29]. This expres¬ 
sion depends on the Reynolds number, density p, velocity u, and 
wire diameter d w , as well as “correlation constants” C f d and Csf 
which are determined empirically. The group used a very coarse 
grid, and a time step that used «100 samples per period of the oscil¬ 
lations (at 53 Hz). Their simulation covered a total time of 1.5 s and 
did not reach a steady state of constant amplitude [27]. 


2.2. Present model of a thermoacoustic engine 

Building upon this previous work, our group has built a model 
of a standing wave, 2/4 resonator thermoacoustic engine. It was 
modeled after an experimental engine that is open to the environ¬ 
ment. Using Gambit, the simulation domain is 150 mm long, 
12 mm high. Thirty millimeters away from the left-hand side, the 
stack channels are located. The stack is 10 mm in length and its 
channels are 0.5 mm high. The mesh is shown in Fig. 2. It contains 
a;l 20,000 triangular cells. This type of cells is chosen in order to be 
able to increase the cell density towards the hydrodynamically rel¬ 
evant stack areas. The stack walls were replicated with the mesh, 
which allows for a more detailed analysis of the flow physics inside 
the stacks during a thermoacoustic cycle than modeling the stack 
using porous formulations and body forces. The area to the left of 
the stacks is referred to as “compliance”. The vertical wall of the 
compliance is the pressure anti-node and the velocity node. 

The horizontal stack walls are defined to have a non-zero thick¬ 
ness and are given a heat transfer coefficient as well as the driving 
temperature gradient. In order to initialize a the flow field in the 
computational domain, the left-hand side of the domain is initially 
defined as a pressure inlet. All walls except the horizontal walls in 
the regenerator are defined as adiabatic. The heat transfer coeffi¬ 
cient on the horizontal regenerator walls resulted in the thermal 
interaction between the walls and the working gas and the replica¬ 
tion of the thermoacoustic effect. To account for turbulence, the k- 
e model is used. This model is chosen because it is the one used 
most widely and is valid for a wide range of flows [30]. The 
description of the model refers to the creation (fe) and dissipation 
(e) of turbulent kinetic energy. Once a steady solution to this con¬ 
dition is found, this area is redefined as a wall and the transient 


Table 1 

Boundary conditions applied to CFD model. 


Parameter 

Value 

Working gas 

Air 

Equation of state 

Ideal gas 

Mean (ambient) pressure 

101,325 Pa 

Initial pressure disturbance 

10 Pa 

Pressure outlet 

0 Pa gauge (ambient) 

Temperature difference (stack) 

Thot - Tcoid = 700-300 K 

Heat transfer coefficient 

50 W/m 2 1< 

Turbulence model 

k-E 

Transient time step 

10 ps 


108000 


106000 
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Fig. 3. Simulation of pressure oscillations due to the thermoacoustic effect a CFD 
analysis as developed by Zink et al. [31], 


simulation is started. All other boundary conditions are left un¬ 
changed. During the transient simulation, the temperature gradi¬ 
ent is the only driving force in the system and is the sole cause 
of the amplification of the pressure amplitude. The absolute pres¬ 
sure at the compliance wall is used as the metric to illustrate the 
onset of strong thermoacoustic oscillations. The applied boundary 
conditions are summarized in Table 1. 

Fig. 3 shows the achieved oscillations through the analysis of 
pressure at the pressure anti-node. Our model is an advancement 
of previous simulations of thermoacoustic engines because it used 
a physical working gas and a heat input as its sole driving force 
(rather than imposed oscillations) [31]. The oscillation frequency 
is 614 Hz, so the transient time step yields a resolution of 164 time 
steps per period. Since the operating frequency is determined by 
the overall length of the device this does not change when the 
cooling stack is added to the model. 



Fig. 2. Simulation domain used for the initial simulation of the thermoacoustic effect by Zink et al. [31], 
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Fig. 4. Simulation domain used for cooling, TAE with a secondary cooling stack. 



Timestep [1e-5 s] 

Fig. 5. Temperature behavior across the cooling stack. 


3. Simulation of thermally driven cooling 

As a continuation of our simulation efforts, we developed a 
model of a thermoacoustic refrigerator. Thermally driven cooling 
requires a secondary stack where the pressure amplitude of the 
sound generated in the driving stack can be attenuated and with¬ 
draw thermal energy from the gas. Again, this phenomenon formed 
the basis of the investigations by Cao and ensuing models as de¬ 
scribed above. In this study, however, this effect is simulated as 
part of a whole engine rather than a single thermoacoustic couple 
that is exposed to an oscillating pressure boundary. 

3.1. Extended CFD model to include a cooling stack 

The present model extends our previous work by the addition of 
a secondary stack just to the right of the driving stack. As in the 
driving stack, the cooling stack is modeled with physical flow 
channels of height 0.5 mm. Its length is 5 mm. The heat transfer 
coefficient in the cooling stack is chosen to be identical to the driv¬ 
ing stack at 50 W/m 2 K. This model design results in a full repre¬ 
sentation of a heat driven thermoacoustic refrigerator that does 
not rely on mechanical excitation through oscillating boundary 
conditions. Furthermore, our model also advances the investiga¬ 
tion of thermoacoustic devices as a whole rather than an individual 
stack plate. In contrast to the previous models that also simulated 
thermoacoustic cooling, our model is capable of achieving this by 
solely applying a heat input to the system. Fig. 4 illustrates the 
new computational domain. The larger driving stack is located clo¬ 
ser to compliance the device, the new refrigeration stack is located 
to its right. 

Recalling that the closed end (on the left-hand side) is the pres¬ 
sure anti-node and the velocity node of the standing wave. The 
thermoacoustic effect requires large pressure variations and non¬ 
zero displacement. The cooling stack is located further away from 



the closed end. As a result, the pressure amplitude encountered 
here is decreased, but the velocity amplitude is higher than within 
the driving stack. These operating conditions are not ideal for the 
operation of a TAR but they can replicate the concept of thermo- 
acoustically achieved cooling. 

3.2. Demonstration of thermally driven cooling 

Fig. 5 shows the resulting temperature behavior throughout the 
cooling stack with respect to time. The variable x* refers to a local 
coordinate system with its origin at the left-hand side of the cool¬ 
ing stack. The thin lines show the detailed gas temperature, the 
thick lines its average value for each period. It is clear that the aver¬ 
age temperature achieved in the cooling stack reaches below the 
ambient temperature of 300 K. Due to the lack of exhausting the 
pumped heat from the hot side of the cooling stack, the effect only 
lasts a short period of time. 

We compare these results with measurements that were taken 
from a thermoacoustically driven refrigerator that was set up as 
shown in the simulation mesh, that is two stacks (utilizing ex¬ 
truded cordierite) located close to each other. The larger stack 
was heated with an electric heating wire. Due to the low thermal 
conductivity of cordierite, active cooling of the cold side of the 
stack is not required. The secondary cooling stack was located close 
to, but not in contact with the driving stack, towards the open end 
of the device. Fig. 6 shows the temperature difference across the 
refrigeration stack. At the onset of oscillations, a temperature dif¬ 
ference between the two stack ends is established. This is due to 
thermoacoustic heat pumping. There was neither a low tempera¬ 
ture load nor an exhaust heat exchanger, which explains why the 
cold side temperature eventually starts rising. The internal heat 
conduction through the cooling stack’s solid areas results in a heat 
flux from the hot side to the cold side. 
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Fig. 7. Velocity distribution throughout the cooling stack, plotted for one entire period. 


3.3. Velocity, temperature, and heat flux with respect to time 

The flow variables of interest in thermoacoustics are the veloc¬ 
ity and temperature, and, as a result, the heat transfer between the 
solid and gas inside the cooling stack. The following data is ex¬ 
tracted from a simulation that reached the state of sustained pres¬ 
sure oscillations, after r; 44,000 time steps. The velocity component 
is shown in Fig. 7. This figure shows the velocity component of the 
wave throughout the cooling stack in the center of the channel 
with respect to x* and time. Interpreting the velocity behavior, 
the undisturbed sinusoidal behavior is obvious. Also, the velocity 
inside the stack is not significantly influenced by the presence of 
the stack. It must be noted that for the period shown, initially, 
the axial velocity is negative; that is, gas is streaming from the res¬ 
onator through the stacks towards the closed end of the device. 
This supports the explanation of the behavior of the temperature 
and the heat flux given below. 

With the (initially negative) velocity component, the interpreta¬ 
tion of the temperature behavior shown in Fig. 8 can be explained. 
As the velocity becomes positive and continues to increase, there is 
a disturbance in the temperature contour that propagates through 
the stack. It dissipates before the opposite end is reached. One half 
period after this disturbance forms on one end of the cooling stack, 
a new disturbance starts on the opposite side. This is caused by the 
reversal of flow and the fluid interacting with the cooling stack. 
This disturbance is the same phenomenon that was discussed pre¬ 
viously by Cao et al. [8] and Marx et al. [20]. The visualization of 
the temperature with respect to time and axial position is a novel 
approach to visualizing the effects of the stack on flow variables 
during the thermoacoustic cycle. 

Next, comparing the temperature behavior of the channel cen¬ 
ter and in the vicinity of the wall, clear difference between the 
two locations become apparent. The aforementioned disturbance 
caused by the entrance of the stack starts later and dissipates soon¬ 
er. Furthermore, the disturbance that develops during the return 
cycle is more developed than in the center of the channel. Of 
course, it is the difference between the two temperatures that is 
proportional to the heat exchange between gas and solid stack 
wall. 




Fig. 8. Gas temperature in the center and on the wall of a channel along the length 
of the cooling stack, plotted for one entire period. 
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Fig. 9. Heat transfer on the wall Q.,.. along the length of the cooling stack plotted for one entire period. A Q= 0 surface is plotted in addition for improved visualization. 


Investigating the heat transfer on the channel wall Qw 
throughout the stack, shown in Fig. 9, the sinusoidal behavior 
with respect to time is apparent. The additional surface cutting 
through the heat transfer surface illustrates Q=0 for improved 
visualization of Q^,. The temperature disturbance discussed above 
affects the heat transfer, which results in a similar disturbance to 
propagate through the surface plot of the heat transfer. In addi¬ 
tion, and more importantly, it is shown that the time averaged 
heat transfer at x* = 5 mm is less than 0, while the time averaged 
heat transfer at the opposite end of the stack at x* = 0 mm is 
greater than 0. A peak is discernible but its visualization is limited 
due to the discretization of the data extraction during the 
simulation. 

4. Summary and outlook 

We have demonstrated the extension of previous modeling ef¬ 
forts of a thermoacoustic engine [31] to include a cooling stack. 
Rather than previous investigations that illustrated thermoacoustic 
cooling with oscillations imposed as a boundary condition, our 
new model only utilizes thermal energy as an input. The oscilla¬ 
tions are thus created within the model as part of the simulation. 
From our previous model, the TAE was driven using a temperature 
gradient imposed along the horizontal stack surfaces. Oscillations 
were induced with a small initial pressure disturbance. Once ther¬ 
moacoustic oscillations were achieved, the temperature of the gas 
was recorded along the cooling stack. A decrease in temperature 
below ambient was proven. This is consistent with temperature 
measurements taken from a similar physical model. Currently, per¬ 
formance is limited by the design of the TAR, specifically the loca¬ 
tion of the cooling stack in the direct proximity of the driving stack 
towards the open end of the modeled 7/4 resonator. This location is 
responsible for relatively high velocities, and thus viscous losses, 
while providing a comparatively small pressure amplitude. Locat¬ 
ing this cooling stack closer to the pressure node should yield bet¬ 
ter performance. Nonetheless, the present model serves as an 
advancement in the demonstration of the principle of thermally 
driven cooling using CFD analysis. 
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